//-------------------------------------------------------------------------//
//                                                                         //
//  This benchmark is an OpenMP C version of the NPB LU code. This OpenMP  //
//  C version is developed by the Center for Manycore Programming at Seoul //
//  National University and derived from the OpenMP Fortran versions in    //
//  "NPB3.3-OMP" developed by NAS.                                         //
//                                                                         //
//  Permission to use, copy, distribute and modify this software for any   //
//  purpose with or without fee is hereby granted. This software is        //
//  provided "as is" without express or implied warranty.                  //
//                                                                         //
//  Information on NPB 3.3, including the technical report, the original   //
//  specifications, source code, results and information on how to submit  //
//  new results, is available at:                                          //
//                                                                         //
//           http://www.nas.nasa.gov/Software/NPB/                         //
//                                                                         //
//  Send comments or suggestions for this OpenMP C version to              //
//  cmp@aces.snu.ac.kr                                                     //
//                                                                         //
//          Center for Manycore Programming                                //
//          School of Computer Science and Engineering                     //
//          Seoul National University                                      //
//          Seoul 151-744, Korea                                           //
//                                                                         //
//          E-mail:  cmp@aces.snu.ac.kr                                    //
//                                                                         //
//-------------------------------------------------------------------------//

//-------------------------------------------------------------------------//
// Authors: Sangmin Seo, Jungwon Kim, Jun Lee, Jeongho Nah, Gangwon Jo,    //
//          and Jaejin Lee                                                 //
//-------------------------------------------------------------------------//

#include <stdio.h>
#ifdef _OPENMP
#include <omp.h>
#endif
#include "applu.incl"
#include "timers.h"

#include <mpi.h>

//---------------------------------------------------------------------
// Thread synchronization for pipeline operation
//---------------------------------------------------------------------
/* common /threadinfo1/ */
int isync[ISIZ2+1];
/* common /threadinfo2/ */
int mthreadnum, iam;
//#pragma omp threadprivate(mthreadnum,iam)


//---------------------------------------------------------------------
// to perform pseudo-time stepping SSOR iterations
// for five nonlinear pde's.
//---------------------------------------------------------------------
void ssor(int niter)
{
  //---------------------------------------------------------------------
  // local variables
  //---------------------------------------------------------------------
  int i, j, k, m, n;
  int istep;
  double tmp, tmp2, tv[ISIZ2][ISIZ1][5];
  double delunm[5];



  //---------------------------------------------------------------------
  // begin pseudo-time stepping iterations
  //---------------------------------------------------------------------
  tmp = 1.0 / ( omega * ( 2.0 - omega ) );

MPI_Status status;
   int my_id = 0, id_num = 0,upper_bound = 0, lower_bound = 0, pros_num = 0, average_span = 0, root_process=0;
   int ierr,upper_bound_two = 0, lower_bound_two = 0, average_span_two = 0;
   int upper_bound_three = 0, lower_bound_three = 0, average_span_three = 0, last_pros_bound_three =0;

///*
int *disps2;
int rcounts2[pros_num];
int *kdeltas2;
 int p2;

int *disps7;
int rcounts7[pros_num];
int *kdeltas7;
int p7;

//*/
   MPI_Comm comm;

   root_process = 0;

   ierr = MPI_Comm_rank(MPI_COMM_WORLD, &my_id);
   ierr = MPI_Comm_size(MPI_COMM_WORLD, &pros_num);
/*
int *disps2;
int rcounts2[pros_num];
int *kdeltas2;
 int p2;
*/
   average_span = (jend-jst)/pros_num;
   average_span_two = (nz)/pros_num;
   average_span_three = ((nz-1)-1)/pros_num;

        if (my_id == root_process){

        lower_bound = jst;
        upper_bound = average_span;

        lower_bound_two = 0;
        upper_bound_two = average_span_two;

        lower_bound_three = 1;
        upper_bound_three = average_span_three;


           if ((my_id+1)  == pros_num){
                if(upper_bound < jend ){
                        upper_bound = jend;
                   }
                if(upper_bound_two < nz ){
                        upper_bound_two = nz;
                   }
                if(upper_bound_three < (nz-2) ){
                        last_pros_bound_three = ((nz-2)-lower_bound_three)+1;
                        upper_bound_three = nz-2;
                   }
                }

         }
        else{

        lower_bound = (average_span *  my_id) + 1;
        upper_bound = average_span * (my_id+1);

        lower_bound_two = (average_span_two *  my_id) + 1;
        upper_bound_two = average_span_two * (my_id+1);

        lower_bound_three = (average_span_three *  my_id) +1;
        upper_bound_three = average_span_three * (my_id+1);


          if ((my_id+1) == pros_num){
                if(upper_bound < jend ){
                        upper_bound = jend;
                   }
                if(upper_bound_two < nz ){
                        upper_bound_two = nz;
                   }
                if(upper_bound_three < (nz-2) ){
                        last_pros_bound_three = ((nz-2)-lower_bound_three)+1;
                        upper_bound_three = (nz-2);
                   }
                 }
             }





  //---------------------------------------------------------------------
  // initialize a,b,c,d to zero (guarantees that page tables have been
  // formed, if applicable on given architecture, before timestepping).
  //---------------------------------------------------------------------
//  #pragma omp parallel default(shared) private(m,n,i,j)
  {
//  #pragma omp for nowait
//  for (j = jst; j < jend; j++) {
 for (j = lower_bound; j < upper_bound; j++) {
    for (i = ist; i < iend; i++) {
      for (n = 0; n < 5; n++) {
        for (m = 0; m < 5; m++) {
          a[j][i][n][m] = 0.0;
          b[j][i][n][m] = 0.0;
          c[j][i][n][m] = 0.0;
          d[j][i][n][m] = 0.0;
        }
      }
    }
  }
//  #pragma omp for nowait
//  for (j = jend - 1; j >= jst; j--) {
  for (j = upper_bound-1; j >= lower_bound; j--) {
    for (i = iend - 1; i >= ist; i--) {
      for (n = 0; n < 5; n++) {
        for (m = 0; m < 5; m++) {
          au[j][i][n][m] = 0.0;
          bu[j][i][n][m] = 0.0;
          cu[j][i][n][m] = 0.0;
          du[j][i][n][m] = 0.0;
        }
      }
    }
  }
  } //end parallel
  for (i = 1; i <= t_last; i++) {
    timer_clear(i);
  }

  //---------------------------------------------------------------------
  // compute the steady-state residuals
  //---------------------------------------------------------------------
  rhs();
  //---------------------------------------------------------------------
  // compute the L2 norms of newton iteration residuals
  //---------------------------------------------------------------------
  l2norm( ISIZ1, ISIZ2, ISIZ3, nx0, ny0, nz0,
          ist, iend, jst, jend, rsd, rsdnm );

  for (i = 1; i <= t_last; i++) {
    timer_clear(i);
  }
  timer_start(1);

  //---------------------------------------------------------------------
  // the timestep loop
  //---------------------------------------------------------------------
  for (istep = 1; istep <= niter; istep++) {
    if ((istep % 20) == 0 || istep == itmax || istep == 1) {
      if (niter > 1) printf(" Time step %4d\n", istep);
    }

    //---------------------------------------------------------------------
    // perform SSOR iteration
    //---------------------------------------------------------------------
//    #pragma omp parallel default(shared) private(i,j,k,m,tmp2) \
                shared(ist,iend,jst,jend,nx,ny,nz,nx0,ny0,omega)
    {
//    #pragma omp master
   if(my_id == root_process){
    if (timeron) timer_start(t_rhs);
     }

    tmp2=dt;


//    #pragma omp for nowait
//    for (k = 1; k < nz - 1; k++) {
    for (k = lower_bound_three; k <= upper_bound_three; k++) {
      for (j = jst; j < jend; j++) {
        for (i = ist; i < iend; i++) {
          for (m = 0; m < 5; m++) {
            rsd[k][j][i][m] = tmp2 * rsd[k][j][i][m];
          }
        }
      }
    }

/*
kdeltas = (int *) malloc(pros_num*sizeof(int));
disps = (int *) malloc(pros_num*sizeof(int));
disps[0] = 0;
rcounts[0] = kdeltas[0]*(ISIZ2/2*2+1)*(ISIZ1/2*2+1)*(5);


for(p=0; p<pros_num; p++) {
rcounts[p] = average_span_three*(ISIZ2/2*2+1)*(ISIZ1/2*2+1)*(5);
 disps[p] = disps[p-1] + rcounts[p-1];

}

*/



kdeltas2 = (int *) malloc(pros_num*sizeof(int));
disps2 = (int *) malloc(pros_num*sizeof(int));
disps2[0] = 0;
//rcounts2[0] = kdeltas2[0]*(ISIZ2/2*2+1)*(ISIZ1/2*2+1)*5;
rcounts2[0] = average_span_three*(ISIZ2/2*2+1)*(ISIZ1/2*2+1)*5;
for(p2=1; p2<pros_num; p2++) {
//rcounts2[p2] = kdeltas2[p2]*(ISIZ2/2*2+1)*(ISIZ1/2*2+1)*(5);
rcounts2[p2] = (average_span_three)*(ISIZ2/2*2+1)*(ISIZ1/2*2+1)*5;
//          if ((pros_num >=4)&& ((my_id+1)== pros_num)){
//          rcounts2[p2] =  last_pros_bound_three*(ISIZ2/2*2+1)*(ISIZ1/2*2+1)*(5);
//                }
disps2[p2] = disps2[p2-1] + rcounts2[p2-1];

}
ierr = MPI_Allgatherv(MPI_IN_PLACE,0,MPI_DATATYPE_NULL,&rsd[1][0][0][0],rcounts2,disps2,MPI_DOUBLE,MPI_COMM_WORLD);

//int scount = average_span_three*(ISIZ2/2*2+1)*(ISIZ1/2*2+1)*(5);
//int scount = 5*(ISIZ2/2*2+1)*(ISIZ1/2*2+1)*(5);
//int rcount = scount;
//..ierr=MPI_Allgather( &rsd[lower_bound_three][0][0][0], scount, MPI_DOUBLE, &rsd[1][0][0][0], rcount, MPI_DOUBLE, MPI_COMM_WORLD);


//..ierr = MPI_Allgatherv(MPI_IN_PLACE,0,MPI_DATATYPE_NULL,&rsd[1][0][0][0],rcounts,disps,MPI_DOUBLE,MPI_COMM_WORLD);

//..int MPI_Barrier( MPI_Comm comm );

//change to allgatherv

/*
    tright = 0;
    if(my_id == 0) {
        MPI_Send(&rsd[upper_bound_three-(average_span_three-1)][0][0][0],scount,MPI_DOUBLE,rankp1,tright,MPI_COMM_WORLD);
    } else {
        MPI_Recv(&rsd[lower_bound_three-average_span_three][0][0][0],rcount,MPI_DOUBLE,rankm1,tright,MPI_COMM_WORLD,&status);
         }

    tleft = 1;
    if(my_id == pros_num-1) {
        MPI_Send(&rsd[lower_bound_three][0][0][0],scount,MPI_DOUBLE,rankm1,tleft,MPI_COMM_WORLD);
    } else {
        MPI_Recv(&rsd[upper_bound_three+1][0][0][0],rcount,MPI_DOUBLE,rankp1,tleft,MPI_COMM_WORLD,&status);
     }

*/


//    #pragma omp master
  if(my_id == root_process){
    if (timeron) timer_stop(t_rhs);
      }
    mthreadnum = 0;
    //mthreadnum = omp_get_num_threads() - 1;
    if (mthreadnum > jend - jst) mthreadnum = jend - jst;
    iam = 0;
    //iam = omp_get_thread_num();
    if (iam <= mthreadnum) isync[iam] = 0;
//    #pragma omp barrier
//      int MPI_Barrier( MPI_Comm comm );

    for (k = 1; k < nz -1; k++) {
      //---------------------------------------------------------------------
      // form the lower triangular part of the jacobian matrix
      //---------------------------------------------------------------------
    
//      #pragma omp master
     if(my_id == root_process){
      if (timeron) timer_start(t_jacld);
      }
      jacld(k);
//      #pragma omp master
    if(my_id == root_process)
      {
      if (timeron) timer_stop(t_jacld);

      //---------------------------------------------------------------------
      // perform the lower triangular solution
      //---------------------------------------------------------------------
      if (timeron) timer_start(t_blts);
      }
      blts( ISIZ1, ISIZ2, ISIZ3,
            nx, ny, nz, k,
            omega,
            rsd, 
            a, b, c, d,
            ist, iend, jst, jend, 
            nx0, ny0 );
//      #pragma omp master
	if(my_id == root_process){
      	   if (timeron) timer_stop(t_blts);
                   }
//        int MPI_Barrier( MPI_Comm comm );

    }
//    #pragma omp barrier
//      int MPI_Barrier( MPI_Comm comm );

    for (k = nz - 2; k > 0; k--) {
      //---------------------------------------------------------------------
      // form the strictly upper triangular part of the jacobian matrix
      //---------------------------------------------------------------------
//      #pragma omp master
     if(my_id == root_process){
        if (timeron) timer_start(t_jacu);
       }
//     int MPI_Barrier( MPI_Comm comm );
      jacu(k);
//      #pragma omp master
     if(my_id == root_process)
      {
      if (timeron) timer_stop(t_jacu);

      //---------------------------------------------------------------------
      // perform the upper triangular solution
      //---------------------------------------------------------------------
      if (timeron) timer_start(t_buts);
      }
 
//     int MPI_Barrier( MPI_Comm comm );
      buts( ISIZ1, ISIZ2, ISIZ3,
            nx, ny, nz, k,
            omega,
            rsd, tv,
            du, au, bu, cu,
            ist, iend, jst, jend,
            nx0, ny0 );
//      #pragma omp master
        if(my_id == root_process){
           if (timeron) timer_stop(t_buts);
                   }
//        int MPI_Barrier( MPI_Comm comm );

    }


//    #pragma omp barrier
//      int MPI_Barrier( MPI_Comm comm );
    //---------------------------------------------------------------------
    // update the variables
    //---------------------------------------------------------------------
//    #pragma omp master
     if(my_id == root_process){
          if (timeron) timer_start(t_add);
          }
//     int MPI_Barrier( MPI_Comm comm );

    tmp2 = tmp;

//    #pragma omp for nowait
//    for (k = 1; k < nz-1; k++) {
   for (k = lower_bound_three; k <= upper_bound_three; k++) {
      for (j = jst; j < jend; j++) {
        for (i = ist; i < iend; i++) {
          for (m = 0; m < 5; m++) {
            u[k][j][i][m] = u[k][j][i][m] + tmp2 * rsd[k][j][i][m];
          }
        }
      }
    }
    } //end parallel

double scount_two = average_span_three*(ISIZ2/2*2+1)*(ISIZ1/2*2+1)*5;
//double scount_two = 4*(ISIZ2/2*2+1)*(ISIZ1/2*2+1)*5;
double rcount_two = scount_two;
//ierr=MPI_Allgather( &u[lower_bound_three][0][0][0], scount_two, MPI_DOUBLE, &u[1][0][0][0], rcount_two, MPI_DOUBLE, MPI_COMM_WORLD);
//ierr=MPI_Allgather( MPI_IN_PLACE, scount_two, MPI_DOUBLE, &u[1][0][0][0], rcount_two, MPI_DOUBLE, MPI_COMM_WORLD);




kdeltas7 = (int *) malloc(pros_num*sizeof(int));
disps7 = (int *) malloc(pros_num*sizeof(int));
disps7[0] = 0;
rcounts7[0] = average_span_three*(ISIZ2/2*2+1)*(ISIZ1/2*2+1)*5;

for(p7=1; p7<pros_num; p7++) {

rcounts7[p7] = average_span_three*(ISIZ2/2*2+1)*(ISIZ1/2*2+1)*5;

//          if ((pros_num >=4)&& ((my_id+1)== pros_num)){
//          rcounts7[p7] =  last_pros_bound_three*(ISIZ2/2*2+1)*(ISIZ1/2*2+1)*(5);
//                }

disps7[p7] = disps7[p7-1] + rcounts7[p7-1];
}
ierr = MPI_Allgatherv(MPI_IN_PLACE,0,MPI_DATATYPE_NULL,&u[1][0][0][0],rcounts7,disps7,MPI_DOUBLE,MPI_COMM_WORLD);




//    int MPI_Barrier( MPI_Comm comm );
/*
int tright;
int tleft;
int rankp1 = my_id+1;
int rankm1 = my_id-1;

    tright = 0;
    if(my_id == 0) {
        MPI_Send(&u[upper_bound_three-(average_span_three-2)][0][0][0],scount_two,MPI_DOUBLE,rankp1,tright,MPI_COMM_WORLD);
    } else {
        MPI_Recv(&u[lower_bound_three-(average_span_three-1)][0][0][0],rcount_two,MPI_DOUBLE,rankm1,tright,MPI_COMM_WORLD,&status);
         }

    tleft = 1;
    if(my_id == pros_num-1) {
        MPI_Send(&u[lower_bound_three][0][0][0],scount_two,MPI_DOUBLE,rankm1,tleft,MPI_COMM_WORLD);
    } else {
        MPI_Recv(&u[upper_bound_three+1][0][0][0],rcount_two,MPI_DOUBLE,rankp1,tleft,MPI_COMM_WORLD,&status);
     }

*/

    if (timeron) timer_stop(t_add);

    //---------------------------------------------------------------------
    // compute the max-norms of newton iteration corrections
    //---------------------------------------------------------------------
    if ( (istep % inorm) == 0 ) {
      if (timeron) timer_start(t_l2norm);
      l2norm( ISIZ1, ISIZ2, ISIZ3, nx0, ny0, nz0,
              ist, iend, jst, jend,
              rsd, delunm );
      if (timeron) timer_stop(t_l2norm);
      /*
      if ( ipr == 1 ) {
        printf(" \n RMS-norm of SSOR-iteration correction "
               "for first pde  = %12.5E\n"
               " RMS-norm of SSOR-iteration correction "
               "for second pde = %12.5E\n"
               " RMS-norm of SSOR-iteration correction "
               "for third pde  = %12.5E\n"
               " RMS-norm of SSOR-iteration correction "
               "for fourth pde = %12.5E\n",
               " RMS-norm of SSOR-iteration correction "
               "for fifth pde  = %12.5E\n", 
               delunm[0], delunm[1], delunm[2], delunm[3], delunm[4]); 
      } else if ( ipr == 2 ) {
        printf("(%5d,%15.6f)\n", istep, delunm[4]);
      }
      */
    }
 
    //---------------------------------------------------------------------
    // compute the steady-state residuals
    //---------------------------------------------------------------------
    rhs();
    //---------------------------------------------------------------------
    // compute the max-norms of newton iteration residuals
    //---------------------------------------------------------------------
    if ( ((istep % inorm ) == 0 ) || ( istep == itmax ) ) {
      if (timeron) timer_start(t_l2norm);
      l2norm( ISIZ1, ISIZ2, ISIZ3, nx0, ny0, nz0,
              ist, iend, jst, jend, rsd, rsdnm );
      if (timeron) timer_stop(t_l2norm);
      /*
      if ( ipr == 1 ) {
        printf(" \n RMS-norm of steady-state residual for "
               "first pde  = %12.5E\n"
               " RMS-norm of steady-state residual for "
               "second pde = %12.5E\n"
               " RMS-norm of steady-state residual for "
               "third pde  = %12.5E\n"
               " RMS-norm of steady-state residual for "
               "fourth pde = %12.5E\n"
               " RMS-norm of steady-state residual for "
               "fifth pde  = %12.5E\n", 
               rsdnm[0], rsdnm[1], rsdnm[2], rsdnm[3], rsdnm[4]);
      }
      */
    }

    //---------------------------------------------------------------------
    // check the newton-iteration residuals against the tolerance levels
    //---------------------------------------------------------------------
    if ( ( rsdnm[0] < tolrsd[0] ) && ( rsdnm[1] < tolrsd[1] ) &&
         ( rsdnm[2] < tolrsd[2] ) && ( rsdnm[3] < tolrsd[3] ) &&
         ( rsdnm[4] < tolrsd[4] ) ) {
      //if (ipr == 1 ) {
      printf(" \n convergence was achieved after %4d pseudo-time steps\n",
          istep);
      //}
      break;
    }
  }

  timer_stop(1);
  maxtime = timer_read(1);
}

